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Abstract Model reduction methods are relevant when the computation time of a 
full convection-diffusion-reaction simulation based on detailed chemical reaction 
mechanisms is too large. In this article, we review a model reduction approach 
based on optimization of trajectories and show its applicability to realistic com- 
bustion models. As most model reduction methods, it identifies points on a slow 
invariant manifold based on time scale separation in the dynamics of the reaction 
system. The numerical approximation of points on the manifold is achieved by 
solving a semi-infinite optimization problem, where the dynamics enter the prob- 
lem as constraints. The proof of existence of a solution for an arbitrarily chosen 
dimension of the reduced model (slow manifold) is extended to the case of real- 
istic combustion models including thermochemistry by considering the properties 
of proper maps. The model reduction approach is finally applied to three models 
based on realistic reaction mechanisms: 1. ozone decomposition as a small test 
case; 2. simplified hydrogen combustion for comparison with another model re- 
duction method; 3. syngas combustion as a test case including all features of a 
detailed combustion mechanism. 
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1 Introduction 

The modeling of chemically reacting flows comprises the interplay between flow 
(convection), diffusion, and chemical reaction processes. This interplay is fairly 
complex if the model is based on a detailed combustion mechanism involving a 
large number of chemically reactive species and reactions. Here complexity reduc- 
tion and model reduction methods can be effective tools. 

A common aim of many model reduction approaches is the identification (and 
computation) of so called slow invariant manifolds (SIM). Many model reduction 
methods are applied to the chemical source term of the system of reaction trans- 
port partial differential equations (PDE), which describe the reactive flow. This 
means, the model reduction method is only regarding a system of ordinary differ- 
ential equations (ODE) modeling the kinetic source term. Trajectories in chemical 
composition space are relaxing to the SIM while converging towards equilibrium. 
In this sense, the SIM represent the slow chemistry for a time scale separation 
between the tangent and normal dynamics. The existence of a SIM is closely re- 
lated to multiple time scales and time scale separation and a spectral gap in the 
eigenvalues of the Jacobian of the chemical source term characterizes the ratio of 
contraction rates in tangent and normal directions. 

In general, certain species will be chosen as "represented" ones for the sim- 
ulation of the reacting flow based on reduced chemistry models; these are also 
called reaction progress variables. These variables parametrize the reduced chem- 
istry model and for these variables the transport PDE are actually solved. The 
values of the remaining unrepresented variables are computed in dependence of 
the represented species by considering a point on the slow manifold parameterized 
by the reaction progress variables. 

Historically, model reduction techniques have been used since the quasi steady 
state assumption (QSSA) and the partial equilibrium assumption (PEA) became 
popular — methods that had to be performed by hand [39]. By contrast, modern 
model reduction methods compute a slow manifold approximation automatically 
without the need of expert knowledge for identification of reactions in partial 
equilibrium and species in steady state within a complex chemical reaction mecha- 
nism. A very popular automatic method is the intrinsic low dimensional manifold 
(ILDM) method, that has originally been published by Maas and Pope in [28] 
in 1992, and its extensions. Also the computational singular perturbation (CSP) 
method by Lam and Goussis [18, 19] is widely applied, e.g. in [30]. Another widely 
applied method, e.g. in [17], is the method of fiamelet generated manifolds [10]. 
For an overview of model reduction methods, see the review [14] and the references 
it contains. 

The scope of this manuscript is a guidance to the application of the optimiza- 
tion based model reduction method as introduced in [20] and future developments. 
The method is applied to chemical combustion mechanisms and results are dis- 
cussed. The outline of this manuscript is as follows. In Section 2, the chemistry 
models regarded here are presented. The optimization problem for identification 
of the SIM is explained in Section 3. A short overview of appropriate numerical 
methods needed for solving the optimization problem is given in Section 4. Results 
of an application to three models are shown and discussed in Section 5. 
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2 Model equations in combustion chemistry 

The model reduction approach discussed here is apphed only to the reaction part 
of the reactive flow model. In this section we review knowledge that can be found 
in text books as e.g. [16,39]. 

2.1 Mass conservation laws 

The general reaction transport equation in the variable ( which can be mass frac- 
tions, temperature, or any variable describing the state of the system depending 
on time t and space x can be written as 

9tC = <s(c) + r(c,9xCa^c), (1) 

where S is the (chemical) source term and T the physical transport operator, i.e. 
convection and diffusion. 

The model comprises rispec chemical species composed by neiem chemical el- 
ements, and the chemical source term 5 obeys the law of elemental mass con- 
servation and an energetic balance which will be regarded in Section 2.2. In the 
following, the variables Zi are given in terms of specific moles, which are defined 
as the amount of species i (n^) divided by the total mass (m) of the system, which 
is the same as the species's mass fractions {wi) divided by its molar mass (Mj): 

_ Ui _ Wj 

m~ Mi' 

In these variables, the mass conservation of each element in the system is 
formulated as 

"spec 

«i = XI ^^^^^ i = 1, . . . , rieiem, (2) 
3 = 1 

where Xij is the atomic composition coefficient — the number of element i in species 
j. There is also a restriction to the choice of the elemental specific moles Zi requiring 
that the mass fractions sum to one 

"clem 

i=l 

where Mj is the molar mass of element i. This is equivalent to the conservation of 
the total mass of the system. 

2.2 Energetic balance 

Energy conservation has to be regarded, too. We consider systems within one of 
the four standard thermodynamic environments, i.e. 

— isothermal and isochoric 

— isothermal and isobaric 

— adiabatic and isochoric (hence isoenergetic) 

— adiabatic and isobaric (hence isenthalpic) 
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systems. Only the temperature is fixed in the isothermal case, while the specific 
enthalpy h 

^spec 

h=^H°iT)zi, (3) 

i=l 

or the specific internal energy e 

e = h-^, (4) 

respectively, are fixed in the adiabatic cases, where we assume an ideal gas mixture, 
and M is the mean molar mass 

M 1 



The molar enthalpy H°{T) of species i is computed by evaluation of so called 
NASA polynomials [7]. 



2.3 Standard systems 

In the optimization problem discussed later in Section 3, the dynamics of the 
system are considered as constraints. Therefore, the reaction ODE system that 
is given by the source term only is discussed in the following. This ODE system 
includes mass action kinetics and incorporates a differential form of the elemental 
mass conservation laws. 

The mass balance in specific moles Zi can be formulated as 

dtZi := ^^Zi = S"^ := i = 1, . . . , rispec- (5) 

In the right hand side of Equation (5), the symbol p refers to the overall mass 
density in the system which is given via 

_ m _ pM 
^~ V ~ RT' 

In the isochoric case, the total mass and volume V are to be known; in the isobaric 
case the total pressure p, the gas constant R, the temperature T, and the mean 
molar mass are necessary. The molar net chemical production rate is denoted by oj 
in Eq. (5). It has to be computed based on a set of chemical elementary reactions 
and their parameters as described in Section 2.4. 

In our case, we formulate the energy balance via the right hand side of the 
temperature equation dtT = = S^. In the isothermal case, we have 5® = of 
course. In the adiabatic cases, energy or enthalpy conservation define the concise 
form of S'^ 
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2.4 Chemical kinetics 

In the remainder of this section, we review the computation of w. A chemical 

combustion mechanism generally is given as a set of rireac elementary reactions 
involving rispec species (and eventually a third body M) 

^spec '^spec 

Z/-jXj ^ ^ t'ijX,, j = 1, . . . , nreac 
i=l i=l 

with the chemical species Xj and the forward and reverse stoichiometric coefficients 
vlj and ul'j . The forward and reverse rate of reaction j is given via 



1=1 

"spei 



i=l 



with the concentrations Ci of species i and the rate constants kfj and k^j . Using 
the net stoichiometric coefficient 



II _ I 

the net rate of change Wj of species i is computed as 

^reac 

Wi = ^ nj{rf,j -rr,j). (7) 
i=i 

In case a third body M takes part in reaction j, third body collision efficiencies aj 
for all species i = 1, . . . , rispec are to be given. The third body concentration 

^spec 

Cm = ^ oiiCi (8) 

i=l 

is then multiplied to the products in Equation (6). 

Formulas for the computation of the rate coefficients are stated in the following. 
The elementary reactions in the mechanism considered here are given in Arrhenius 
format and pressure dependent Troe format, resp. The three parameters A, b, and 
-Ea are given in the Arrhenius kinetics for each reaction. The forward reaction rate 
coefficient is computed via the extended Arrhenius formula 

kij = ^ TiK e~RT. (9) 

A more complicated formula applies in case of pressure dependent reactions. 
Here k{j is computed using Troe fall off curves [12,37]. Both the low pressure 
rate coefficient fco and the high pressure coefficient koo are given via the extended 
Arrhenius formula (9). These are used together with a third body M to compute 
the reduced pressure 

fco Cm 

= -r. — ' 
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where cm is defined as in Equation (8). With this parameter, the final rate constant 
is computed as 



l+Pr 



with a function F. To compute the value of F, Gilbert, Luther, and Troe [12,37] 
introduced the formula 



lgF= <j 1 + 
with a set of simplifications 



IgPr + C 



n - d(lgpr + c) 



c = -0.4 -0.67 IgFc 
n = 0.75- 1.271gFc 
d = O.U 



and the F-center-value 

Fc = (1 - a) exp 



T 



+ a exp 



T 



+ exp 



T* 



which includes the four parameters a, T* , T**, and T***. These are given for each 
so called Troe-reaction. 

The reverse reaction rate constant kfj of reaction j is computed via the equi- 
librium constant K^j of the reaction.^ The equilibrium constant of reaction j in 
terms of concentrations is given as 



P 

RT 



exp 



RT 



with the standard pressure p° and the net change of the number of species present 
in the gas phase 



The change in entropy AS^j and enthalpy AH^-j can be computed by using an 
evaluation of the NASA polynomials for their molar values of the species involved. 
The reverse rate of reaction j finally is 



^ In some publications, the reverse rate coefficients of Arrfienius type reactions are com- 
puted witfi fitted Arrfienius parameters for the reverse reactions. Wc also used this strategy 
in previous publications as e.g. [22,23]. However, this is impossible in case of pressure depen- 
dent reactions and it may lead to inconsistent values of thermodynamic quantities in the mass 
action kinetics in relation to the heat of reaction. Therefore, we prefer the thermodynamic 
approach in this manuscript. 
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3 Optimization problem 

In 2004, Lebiedz introduced a model reduction method based on the minimization 
of entropy production rate along trajectories in chemical composition space [20]. 
The basic idea is that the SIM is characterized by maximum relaxation of the 
system dynamics under given constraints of fixed reaction progress variables. This 
approach has been extended and refined in a number of following publications [21, 
22,25,33]. 

The general geometric situation in the phase space of the reaction system 
spanned by the variables can be seen in the sketch in Figure 1. Trajectories 
bundle on the manifolds of slow motion, that are hierarchically ordered. The aim 




Z2 



Fig. 1 Sketch of the chemical composition space. The domain, where the dynamics take 
place is the blue-bounded polytope, here three-dimensional. Within this polytope there might 
be (depending on a possible time scale separation) a two-dimensional manifold (depicted in 
red), where (shown in green) trajectories relax onto. The trajectories relax then onto the one- 
dimensional manifold within the two-dimensional one. Finally the trajectories converge toward 
the zero-dimensional manifold: the equilibrium. 



is to compute an approximation of such a manifold of given dimension point-wise 
such that the free variables are computed depending on the given reaction progress 
variables which parametrize the manifold. 

Following the idea of [20-22,25,33], an optimization problem has to be solved 
for the approximation of points on the manifold. It can be written in specific moles 
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and temperature as minimization of an objective functional 



min / ^{z{t)) 



(10a) 



subject to 



dtT{i) 



= S-(^(t),T(t)) 
^) = S%z{t),T{t)) 

= C{z{u),T{u)) 



(10b) 
(10c) 

(lOd) 
(lOe) 
(lOf) 







= 2:j(t*) - 0**, jell 
Q^z{t),T{t) 



-rpv 



and 



t € [toM 

U€[to,tf] (fixed). 



(lOg) 
(lOh) 



In the following, we explain optimization problem (10) in detail starting with the 
constraints. 

Chemical source and heat of reaction For the numerical solution of the optimization 

problem, the dynamics of the system have to be computed. The dynamics are given 
via the ODE system in the constraints (10b) and (10c). This ensures to identify 
as solution of (10) a special solution trajectory piece. 

Conservation relations All necessary additional conservation laws are combined in 
the (nonlinear) function C. As discussed before, the dynamics (10b) and (10c) 
contain differential forms of the balances of mass and energy. The concise values 
of the conserved quantities have to be specified at some (fixed) point in time along 
a solution which we choose to be t*. This is Equation (2) and a specification of a 
fixed temperature, enthalpy, or energy, as e.g. (3) or (4), depending on the assumed 
thermodynamic environment. 

SIM parametrization In order to approximate the SIM, a parametrization needs 
to be specified. The species (i.e. the specific moles Zj) which serve as reaction 
progress variables and especially their number have to be specified in advance. 
The number of progress variables determines the chosen dimension of the SIM to 
be approximated. 

In the case illustrated in Figure 1, one might choose zi and Z2 as reaction 
progress variables for parametrization in order to compute a value for the remain- 
ing free variable 2:3 which is supposed to be on the two-dimensional manifold. 
Alternatively, one can choose only zi as reaction progress variable for approxima- 
tion of the one-dimensional manifold. 

The indices of the reaction progress variables are collected in the index set 
Irpv C {1, . . . , rispec}, and their values are fixed in the optimization problem to 
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Positivity Since only positive values of specific moles and temperature have a phys- 
ical meaning, this is included in (lOf). In the optimization context with a realistic 
combustion mechanism included in the constraints, this is also technically impor- 
tant as negative values of Zi and T can result in undefined values (logarithm of a 
negative number) of the right hand sides S"" and 5"^ . The positivity and the linear 
mass conservation relations define a polytope as depicted in Figure 1. 

The chosen point in time U £ [to , if] specifies the position where the fixation of 
the reaction progress variables and the constraint C is applied along the trajectory 
piece, which is optimized. In first publications, e.g. [20,22], = to is chosen. This 
incorporates the demand that the trajectories are fully relaxed to the SIM at time 
t* and no further relaxation takes place afterward. 

The inverse idea is the fixation at the end point t» = tf. The solution point 
and solution trajectory piece of the optimization problem (10) is supposed to be 
part of the SIM. Therefore, also in backward direction of time the trajectory is 
supposed to be already relaxed. 

This is related to the definition of positive and negative invariance of a set 
under a flow in dynamical systems theory. 

Definition 1 (Invariant set, [40]) Let Q C R" be a set. It is called invariant under 
the vector field C = 5(C), C e if for any Co € ^2 it holds that C(i;Co) G ^ for 
all t e K, where C(i; Co) denotes the solution of the initial value problem C = 5(C) 
with initial value Co at t = 0. It is called positively invariant if this conditions holds 
for positive t ^ and negatively invariant if the conditions holds for negative t ^ 0. 

Clearly, the essential degrees of freedom of the optimization problem is the ef- 
fective phase space dimension rispec — neiem minus the number of reaction progress 
variables |Irpv|- The goal of solving the optimization problem (10) is the deter- 
mination (species reconstruction) of the "missing" values Zi{t*), i ^ Irpv, as a 
function of the parameters z^*, j € Trpv 



3.1 The objective functional 

The relaxation criterion <f is supposed to measure the degree of chemical force 
relaxation along a trajectory. Several criteria have been tested for their SIM ap- 
proximation quality, especially in [22]. 

The SIM to be approximated is considered to be slow. This means, the residence 
time of the trajectory in some open neighborhood of a point on the SIM should 
be large — conversely the change and (to second order) the rate of change of the 
variable values is supposed to be small. A similar idea has been pointed out already 
in [13]. The rate of change is closely related to the curvature of the trajectories 
as geometrical objects in phase space. The rate of change of the specific moles is 
simply z = 5^". Its rate of change is the second derivative 

m = Js-^izitim) s^iz{t),T{t)), 

where we denote the Jacobian of a function S a.s Jg. This can be seen as a direc- 
tional derivative of the chemical source w.r.t. its own direction ii which should be 
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regarded normalized v := = y^^y^ 

where || • ||2 denotes the Euchdean norm. The evaluation of this expression within 
the integral should be done in arc length. A re-parametrization cancels out the 

norm \\S^{z{t),T{t))\\2 such that (in notation that coincides with the general prob- 
lem (10)) a reasonable candidate for the criterion $ would be 

Hz{t)) = \\Js^izit),Tit)) 5-(z(t),r{t))|il. (11) 



3.2 Solution of the optimization problem 

In [25] , the authors study theoretical properties of the optimization based model 
reduction method as described in the sections before. It is shown there always exists 
a solution of the optimization problem (10) with only linear (mass conservation) 
constraints if there exists a feasible solution. 

In case of a realistic combustion mechanism as a model for an adiabatic system, 
the nonlinear internal energy conservation or enthalpy conservation comes into 
play. The existence proof will be extended to this cases in the following. The 
crucial point is the compactness of the feasible domain, which is more complicated 
to ensure in the nonlinear case. 

A simple way to guarantee compactness of the feasible domain would be an 
upper bound for the temperature. Together with the compactness argument for 
the linear constraints [25], the compactness of the feasible domain is obvious. But 
it is not clear at all where to choose the upper cut off for the temperature. 

We avoid this temperature cut off but make use of the definition of molar 
enthalpy via NASA polynomials. Thereby we accept the temperature to outrange 
the domain where the NASA polynomials approximate the molar enthalpy of the 
species appropriately. 

The specific enthalpy /i of a system is given as 

"spec 

h=J2HHT)zi. (12) 

i=l 

The equation for the specific internal energy e is 

"spec 

e = h-RTJ2z^■ (13) 

i=l 

The molar enthalpy H° (T) of species i is a continuous function in the temperature 
T. In our case, it is computed by evaluation of the NASA polynomials. Their 
formula is 

Him =ae + a,T + «f + ^^T^ + "^T^ + '^T^ (14) 
-ft z o 4 o 

with two sets of coefficients a^, i = 1,...,6. One set is given for a temperature 

lower than a certain switch temperature T < Tsw and one set of coefficients for high 

temperature T > Tsw The two branches are connected at Tsw at least continuously. 

There are also upper and lower bounds for the temperature, where the polynomial 

approximation is valid. We ignore these bounds for the following theory. 
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Definition 2 (Proper map, [26]) Let X and Y be topological spaces. A map 
(continuous or not) H : X ^ Y is called proper if the preimage H~^{K) of each 
compact subset K cY is compact. 

To formulate a sufficient condition for properness we need the 

Definition 3 (Divergence to infinity, [26]) If X is a topological space, a sequence 
(xn) in X is said to diverge to infinity if for every compact set if C X there are at 
most finitely many indices n with element Xn € K. 

A sufficient condition for properness is the following 

Lemma 1 (Properness condition, [26]) Suppose X and Y are topological spaces, 
and H : X Y is a continuous m,ap. If X is a second countable Hausdorff space and 
F takes sequences diverging to infinity in X to sequences diverging to infimty in Y, 
then F is proper. 

Proof See [26, p. 119]. 

Lemma 2 (Properness of h and e) The specific enthalpy h and the specific internal 
energy e defined via NASA polynomials seen as functions in T and z are proper maps. 

Proof The vector space K"sp»<:+i is a second countable Hausdorff space. Any non- 
constant polynomial takes sequences diverging to infinity in M" equipped with its 
Euclidean metric induced topology to sequences diverging to infinity in R. We can 
see h and e as polynomials of sixth degree in Zi and T, see Eq. (12), (13), and (14). 
Therefore, h : M"=p'"=+i K and e : M^^p^^+i -j. M are proper maps. □ 

Using this information, we can extend the existence lemma 2.1 in [25]. 

Lemma 3 The feasible set at 

M = {{z, T) : C{z, T) = 0; Zj - zf =Q,je Trpv; {z, T) ^ 0} 

is compact. 

Proof Case 1: isothermal combustion 

The mass conservation together with the positivity and the fixed temperature 
define a polytope in K"spec+i which is closed and bounded, hence (Heine-Borel 

theorem) compact. 

Case 2: adiabatic combustion 

As in the isothermal case, the variables Zi are restricted to a compact polytope 
due to elemental mass conservation and positivity constraints. Following Lemma 2, 
the preimage of the singleton of the fixed energy /enthalpy is a compact subset of 
]]j«spec + i This subset may only be further constrained by the polytope defined by 
the mass conservation and positivity, and the intersection of compact subsets is 
compact. n 

Lemma 4 (Existence of a solution) If the map $ : M"=p««' W in the objective 
functional of the optimization problem (10) is a continuous function and the feasible 
set is not empty, there exists a solution of problem (10). 

Proof Following the argumentation in [25], the semi-infinite optimization problem 
(10) can be reduced to a finite dimensional optimization problem by construction 

of a continuous map {z,T){t^) n> f^' <I> {z(t)) dt. As seen in Lemma 3, the feasible 
set M is compact. Therefore, existence follows from the Weierstrafi theorem. □ 
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4 Numerical methods 

The semi-infinite optimization problem (10) can be solved after suitable discretiza- 
tion of the ODE constraints e.g. either by a sequential quadratic programming 
(SQP) [32] or an Interior Point (IP) method, see e.g. the review [11]. 

4.1 Discretization 

In general, there are two ways for discretization and solution of (10): the sequential 
and the simultaneous approach. 

4.1.1 Sequential approach 

In the sequential approach, ODE solution and optimization are fully decoupled. 
The initial values {z{to),T{to)) are used as optimization variables. Starting at to, 
the system is integrated with a stiff ODE solver, e.g. via a backward differentia- 
tion formulae (BDF) scheme [8]. The integrand in the objective function (10a) is 
integrated itself, and the end point is evaluated in sense of a Mayer term objective 
functional. The optimization iteration is performed after that based on the results 
of the integration and computed derivative information. In the cases U = to and 
t* = tf, a single shooting is appropriate as we deal with a stable ODE system; 
whereas in case of t* € (to,tf), a double shooting is needed for the values at t*. 

4.1.2 Simultaneous approach 

Sometimes (e.g. in case of unstable or extremely stiff systems), it is beneficial 
to use an all-at-once approach using collocation formulae [4]. In this simultane- 
ous approach, the solution of the dynamic constraints and the optimization are 
coupled. The interval [to, if] is divided into sub-intervals. Via e.g. a collocation 
method, polynomials are constructed on each sub-interval tangent to the vector 
field of the dynamics (10b) and (10c) approximating their solution, and the cor- 
responding formulae are treated as constraints in the optimization iteration. In a 
collocation approach, we use a Gaufi-Radau formula with linear, quadratic, and 
cubic polynomials, respectively, because they have stiff decay [4]. 

4.2 Solution of the finite-dimensional optimization problem 

In both cases (sequential and simultaneous), the result of the discretization is 
a finite-dimensional nonlinear programming (NLP) problem. This can be solved 
using an SQP algorithm or an IP method. The SQP algorithm treats the inequality 
constraints using an active set strategy, see e.g. [31]. Newton's method is applied 
to the first order optimality conditions of a quadratic approximation of the NLP 
problem including only equality and active inequality constraints. By activating 
and deactivating constraints, the active set in the solution is identified. In contrast 
in an IP method, the inequality constraints are coupled to the objective function 
via a barrier term forcing the iterates into the interior of the feasible domain. The 
resulting equality constrained NLP problem is solved with a homotopy method: 
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Table 1 Ozone decomposition mechanism for the forward rates as in [29]. CoUision efficiencies 
in reactions including M: ao = 1-14, Q02 = 0.40, aos = 0.92. 



Reaction 




A 1 cm, mol, s 


h 


Ea / kjmol-i 


O + O + M 


^ O2 + M 


2.9 X IQi'^ 


-1.0 


0.0 


O3 + M 


^ O + O2+M 


9.5 X 10" 


0.0 


95.0 


O + O3 


^ O2 + O2 


5.2 X 10^2 


0.0 


17.4 



Newton's method is applied to the first order optimahty conditions. In the progress 
of optimization, the barrier parameter is driven to zero to follow an homotopy path 
to the solution of the NLP problem. 



4.3 Algorithms and software 

In Section 5, we present results of an application of the model reduction method. 
We use IPOPT [38] as the main optimization tool. It turned out to be a robust IP 
algorithm appropriate for our problems. For the solution of linear equation systems 
within the optimization algorithm, HSL routines [15] and MUMPS [3], resp., are 
used. Derivatives needed for the optimization are computed with the open source 
automatic diflFerentiation package CppAD [5,6]. 

For discretization of the optimization problem, we use a collocation approach 
based on a Gaufi-Radau method [4]. Alternatively, we use a shooting approach 
including a BDF integrator that has been developed by D. Skanda for [35]. For 
the numerical solution strategies, see also [24]. We use MATLAB for plotting. 



5 Results 

In this section, results for the application of the optimization based model reduc- 
tion method are shown. As this manuscript is focused on the application to realistic 
systems, we skip a discussion of test equations for model reduction methods. These 
can be found in [22,25]. We only consider combustion mechanisms providing the 
complete kinetic and thermodynamic data. 

In [25], it has been shown that the reverse mode (t* = ti) of the method 
identifies the correct SIM in case of a linear test model and the Davis-Skodje test 
model [9], which has an analytically given one-dimensional SIM, for infinite time 
horizon Iq —00. Hence we use the reverse mode for all results presented here. 



5.1 Ozone decomposition 

As a first small test problem, we consider an ozone decomposition mechanism 
including only three allotropes of oxygen, namely atomic oxygen, dioxygen and 
ozone. The mechanism is given in Table 1. 

The thermodynamic data is used in form of NASA polynomial coefficients. 
In comparison to the results in [22], we set up the mechanism in the framework 
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Fig. 2 Results (onc-dimcnsional SIM) for the ozone decomposition mechanism modeled within 
an isothermal and isochoric environment. The value of zq serves as reaction progress variable. 
The optimization problem is solved several time for different values of Zq varying between 
zero and the largest possible value zq- We use the reverse mode (4f = t») with an integration 
interval of tf — to = 10"'' s. The free Zi (in mol/kg) and temperature are plotted versus zq. 



described in Section 2. Reverse rate coefficients are derived from equifibrium ther- 
modynamics. For mass conservation, the elemental specific mole has to be given 
which is 



1000 



15.999 gmor 



62.5 mol kg 



We consider a density of p = 0.2 kg m ^ in the isochoric case and a pressure of 
p = 10'' Pa for isobaric conditions, respectively. 

In case of the ozone mechanism (Table 1), it is technically not necessary to 
demand positiveness of specific moles and temperature because no pressure de- 
pendent reactions are present. Therefore, the mechanism can be evaluated also in 
nonphysical regions (negative species concentrations visible in the plots). 

Results for the four different thermodynamic environments are shown in Fig- 
ure 2, 3, 4, and 5. The model has two degrees of freedom; we compute a numerical 
approximation of a one-dimensional SIM. The (blue) open rings in the plots are 
the solution points {z,T){u). Orbits through these points in forward and reverse 
direction are also shown (blue curves), where the reverse part coincides with the 
optimal trajectory piece {z,T){t), t € [to, if]- The trajectories converge toward 
equilibrium which is shown as full (red) dot on the left hand side of the subfigures. 

In all plots (Figure 2-5), it can be seen that near equilibrium very good results 
can be achieved as the SIM approximation is nearly invariant, i.e. all open dots 
are lying along one (slow) trajectory. But far from equilibrium at large values of 
the reaction progress variable zq, the full dynamics are active, so the invariance 
of the SIM is poor due to a lack of time scale separation. A short relaxation phase 
can be stated, but at least the values are in a reasonable range. 
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Fig. 3 Results (one-dimensional SIM) for the ozone decomposition mechanism modeled within 
an isothermal and isobaric environment. The plot is arranged as in Figure 2. Again we use the 
reverse mode (tf = t*) with an integration interval of tf — to = 10~^ s. 



6 
5 

04 



40 




50 60 

Zo 



70 



2 
1.8 
1.6 
1.4 

CO 

01-2 

K) 

1 

0.8 
0.6 



xlO 



0.4^ 
40 



50 60 

Zo 



5000 
4500 
4000 
3500 



3000 



2500 
2000 



70 



1500^ 
40 




50 60 

Zo 



70 



Fig. 4 Results (one-dimensional SIM) for the ozone decomposition mechanism modeled within 
an adiabatic and isochoric environment. The plot is arranged as in Figure 2. As in the adiabatic 
case, the reaction evolves faster the integration horizon is reduced to t( — tg = 10~* s. 



5.2 Simplified hydrogen combustion mechanism 



In this section, we review a test case for a simplified combustion mechanism. The 
presentation and results are similar to those presented in [25]. 

The reaction mechanism is given in Table 2. In [25], we show a comparison 
to the results of Al-Khateeb et al. in [1,2]. Hence we use the thermodynamical 
data (in form of NASA coefficients) we received from J. M. Powers and A. N. Al- 
Khateeb. The mechanism itself has been published originally in [27]. The simplified 
version shown in Table 2 has been used by Ren et al. in [34]. The mechanism itself 
consists of five reactive species and inert nitrogen, where in comparison to a full 
hydrogen combustion mechanism the species O2, HO2, and H2O2 are removed. 
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Fig. 5 Results (one-dimensional SIM) for the ozone decomposition mechanism modeled within 
an adiabatic and isobaric environment. The plot is arranged as in Figure 2. Again the integra- 
tion horizon is tf — to = 10~** s. 



Table 2 The simplified mechanism as used in [34]. Collision efficiencies M: ajj = l-0iCm2 = 
2.5, aoH = 1.0, «o = 1.0,OH2O = 12.0, 0N2 = l-O- 



Reaction 






A 1 (cm, mol, s) 


h 


/ kJmol-i 


O-hHa 




H + OH 


5.08 X 10"'' 


2.7 


26.317 


H2-HOH 




HaO-t-H 


2.16 X 10°* 


1.5 


14.351 


+ H2O 




2 OH 


2.97 X 10°'^ 


2.0 


56.066 


H2 -l-M 




2H-I-M 


4.58 X 10^'^ 


-1.4 


436.726 


0-l-H + M 




OH -l-M 


4.71 X 10^** 


-1.0 


0.000 


H -1- OH -1- M 




H2O + M 


3.80 X 10^2 


-2.0 


0.000 



The species are involved in six Arrhenius type reactions, where three combina- 
tion/decomposition reactions require a third body for an effective collision. 

Al-Khateeb et al. identified a one-dimensional SIM for a model including this 
mechanism [1]. The model additionally involves the following parameters: The 
combustion is considered in an isothermal and isobaric environment with a tem- 
perature of T = 3000 K and a pressure of p = 101 325 Pa. The elemental mass 
conservation is given in terms of amount of species; it is 

nn -h 2 + "n-OH + 2 riHjO = 1.25 x 10^^ mol 
JT-OH + ?io + "■H2O = 4.15 X 10~^ mol 
2nN2 = 6.64 X 10"^ mol. 

Therefore, the total mass in the system is m = 1.01 x 10~* kg. We continue to use 
the specific moles as our standard variables here and use = ^ in the following. 

The results are shown in Figure 6. There is a very good agreement of our 
results with theirs. Even on both sides of equilibrium, our approximations coincide 
with the correct one-dimensional SIM on its both branches which consist of two 
heteroclinic orbits in this case. 
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^1 

Fig. 6 Three-dimensional plot of the results for the numerical approximation of a one- 
dimensional SIM of the simplified combustion mechanism computed with the reverse mode, 
2H2O reaction progress variable, and tf — ig = lO^'^s. The illustration is similar to Figure 
9 in [1]. The blue bounded polytope shows the physically feasible state space. Green curves 
correspond to some "arbitrary" trajectories for illustration. The red curve depicts the two 
branches of the SIM as derived in [1]. The open red dot represents the unstable fixed point 
{-Re in [1]); the full red dot represents the equilibrium (-R7). Our results are included as blue 
crosses. The state of the species is given as z\ = zjjj i ^2 = ^Oi S'lid ^3 = -^1120 ™ molkg"'^. 

5.3 Syngas combustion mechanism 

As a last example of a full detailed chemistry combustion mechanism, we use a 
syngas combustion extracted from the GRI 3.0 mechanism [36]. It consists of ah 
33 reactions of the GRI 3.0 mechanism which involve no other species than O, O2, 
H, OH, H2, HO2, H2O2, H2O, N2, CO, and CO2. Those 33 reactions can be split 
up into 31 Arrhenius-type and two pressure-dependent ones. 
The overall reaction can be stated as 

H2 + CO + 02^H20-FC02, 

where N2 only serves as a collision partner. We assume a stoichiometric mixture 
of syngas with air in an adiabatic and isochoric environment. As fixed mass den- 



18 



Dirk Lebiedz, Jochen Siehr 




2 C02 4 2 



H20 



O 

Kl 1.5 



2 H20 



■ C02 



Z H20 



* C02 




ZC02 U ^H20 



2H20 * , 




H20 




ZC02 4 2 



H20 



z 

N 22 



2 H20 2 



■ C02 



3000 
2500. 



^ 2000 
1500, 



2H20 ,J „ 



2H20 2V ^ 



C02 



Fig. 7 Visualization of tlie results (two-dimensional SIM) of the model reduction method 
applied to the syngas combustion model described in Section 5.3. The reverse mode is applied 
with a time horizon of lO"'^ s. We approximate a two-dimensional manifold and use the overall 
products 2H9 S'lid 2:002 ^ reaction progress variables. 



sity, we use p = 0.3kgm~^. We assume a ratio of nn^ : nco = 1 ^ Ij and a 
ratio of '■ ~ 1 ■ 3.76. This leads to a unburned mixture of zqq = = 
5.973molkg^^ and z-f^^ = 22.46molkg^^. The specific internal energy of this mix- 
ture at a temperature of T = 1000 K is used as a fixed specific internal energy for 
SIM computation. 

Results for this model are shown in Figure 7. The same style as in Figure 3 
is used. The resulting SIM approximation points (solutions of the optimization 
problem (10)) are shown as open blue dots together with trajectories emanating 
from these and converging to equilibrium, which is shown as full red dot. In the 
two-dimensional case, invariance cannot be seen by eye inspection, but a reasonable 
manifold is found. 



6 Conclusions 

An optimization method is presented that allows for efficient model reduction of 
realistic combustion models. It is applied to three models for testing its applicabil- 
ity. It can be seen that for an appropriate time scale separation the solution of an 
optimization problem approximates points on a slow invariant manifold. The ap- 
plication of the model reduction approach to a realistic syngas combustion model 
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considered in an adiabatic and isochoric environment demonstrates the applica- 
bility to realistic large scale mechanisms. 

Further research will be needed for an identification of an appropriate number 
and choice of the reaction progress variables for large scale mechanisms. 
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